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Abstract. We have developed a new simulation code aimed at studying the stellar dynamics of a galactic central 
star cluster surrounding a massive black hole. In order to include all the relevant physical ingredients (2-body 
relaxation, stellar mas s spectrum, co lli sions, tidal di sruption, ... ), we chose to revive a numerical scheme pioneered 
by Henon in the 70's (Henon 1971b, ^; Henon 1973). It is basically a Monte Carlo resolution of the Fokker-Planck 
equation. It can cope with any stellar mass spectrum or velocity distribution. Being a particle-based method, it 
also allows one to take stellar collisions into account in a very realistic way. This first paper covers the basic version 
of our code which treats the relaxation-driven evolution of a stellar cluster without a central BH. A technical 
description of the code is presented, as well as the results of test computations. Thanks to the use of a binary tree 
to store potential and rank information and of variable time steps, cluster models with up to 2 x 10^ particles can 
be simulated on a standard personal computer and the CPU time required scales as A'p In(A'p) with the particle 
number A'^p. Furthermore, the number of simulated stars needs not be equal to A'p but can be arbitrarily larger. 
A companion paper will treat further physical elements, mostly relevant to galactic nuclei. 
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1. Introduction 

The study of the dynamics of galactic nuclei is likely to be- 
come a field to which more and more interest will be dedi- 
cated in the near future. Numerous observational studies, 
using high resolution photometric and spectroscopic data, 
point to the presence of "massive dark objects" in the cen- 



tre of most bright galaxies (see reviews by 1 


■Cormendy & 


Richstone 1995; Rees 1997; Ford et al. 199f 


i; I 


Vlagorrian 


ct al. 1998; Ho 1999; van der Marel 199£; de Zeeuw 2001), 


including our own Milky Way (Eckart & Gcnzcl 1996, 


1997; Genzel et al. 1997; Ghez et al. 199? , 19991 


,a; Genzel 



et al. 2000 ) . As the precision of the measurements and the 
rigor of their statistical analysis both increase, the con- 
clusion that these central mass concentrations are in fact 
black holes (BHs) with masses ranging from 10^ Mq to 
a few 10^ Mq becomes difficult to e vade, at least in the 
most conspicuous cases ( Maoz 1998| ). Hence, the time is 
ripe for examining in more detail the mutual influence be- 
tween the central BH and the stellar nucleus in which it 
is engulfed. 
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A few particular questions we wish to answer are the 
following: 

— What is the stellar density and velocity profile near 
the centre? 

— What are the stellar collision and tidal disruption 
rates? Were they numerous enough in past epochs to 
contribute significantly to the growth of the central 
BH? Do they occur frequently enough at present to be 
detected in large surveys and thus reveal the presence 
of the BH? 

— What dynamical role do these processes play? Are they 
just by-products of the cluster's evolution or are there 
any feed-back mechanisms? 

— What is the long-term evolution of such a stellar sys- 
tem? Will it experience core collapse, re-expansion, or 
even gravothermal oscillations like globular clusters? 

To summarize, our central goal is to simulate the evo- 
lution of a galactic nucleus over a few 10^ years while tak- 
ing into account a variety of physical processes that are 
thought to contribute significantly to this evolution or are 
deemed to be of interest for themselves. Here is a list of the 
ingredients we want to incorporate into our simulations: 

— Relaxation induced by 2-body gravitational encoun- 
ters. The accumulation of these weak encounters leads 
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to a redistribution of stellar orbital energy and angu- 
lar momentum. The role of relaxation has been exten- 
sively studied in the realm of globular clusters where 
it leads to the slow built-up of a diffuse halo and an in- 
crease in the central density until the so-called "gravo- 
thermal" catastrophe sets in (see Sect. 1.1). Of partic- 
ular relevance to galactic nuclei are early studies that 
demonstrated that, when a central star-destroying BH 
is present, relaxation will lead to the formation of a 



qua 
I97E 


d-stationary cusp in the stellar density ( 


Peebles 


\ Shapiro & Lightman 1976]; 


3ahcaU & Wolf 1976; 


Dokuchaev & Ozernoi 1977a It|; 


_jiglitman & Shapiro 


197' 
197* 


I, Shapiro & Marchant 197 




Cohn & Kulsrud 





Not only can this be a significant source of fuel for 
the central BH, but it can also play a role in the 
cluster's evolution, as specific orbits are systematically 
depleted. The destruction of deeply bound stars may 



lead to heating the stellar cluster ( Marchant fc Shapiro 
198§). 

Stellar collisions. As the stellar density increases in the 
central regions, collisions become more and more fre- 
quent. By comparing the relaxation time and the col- 
lision time (Binney & Tremaine 1987, for instance), 



one sees that collisions can catch up with relaxation 
when Uy > (lnA)^^*14 w 1000 kms~^, where is 
the velocity dispersion in the stellar cluster. In A the 
Coulomb logarithm and K the escape velocity from 
the surface of a typical main sequence (MS) star. But 
even in colder clusters, collisions may be of interest 
as a channel to create peculiar stars (blue stragglers, 
stripped giants, ... ). 

Stellar evolution. As stars change their masses and 
radii, the way they are affected by relaxation, tidal 
disruptions and collisions is strongly modified. Of par- 
ticular importance is the huge increase in cross sec- 
tion during the giant phase and its nearly vanishing 
when the star finally evolves to a compact remnant. As 
they are not prone to being disrupted, these compact 
stars may segregate towards the centre-most part of 



the n ucleus where they may dominat e the density ( Lee 
199|; iMiralda-Escudc fc Gould 200C| , amongst others). 
Furthermore, the evolutionary mass loss (winds, plan- 
etary nebulae, SNe ... ) may strongly dominate the 



BH's g rowth if this gas sinks to the centre (David et al 



198|7b|; [Norman fc Scovillc 1988| ; [Murphy et al. 1991]) . 
BH growth. Any increase of the BH's mass may ob- 
viously lead to a higher rate of tidal disruptions. It 
also imposes higher stellar velocities and, thus, in- 
creases the relative importance of collisions in compari- 
son with relaxation. A further contributing effect is the 
central density built-up due to the adiabatic adapta- 



tion of stellar orbits to the deepening potential ( 


Young 


198|0|; 


Lee & Goodman 198£ ; 


CipoUina fc Berlin 1994; 


iQuinlan et al. 1995|; 


Sigurdsson et al. 1995; 


Leeuwin fc 



Although we restrict our discussion to these physical 
ingredients in the first phase of our work, there are many 
others of potential interest; some of them are mentioned 
in Sect. |^ 

Thus, the evolution of galactic nuclei is a complex 
problem. As it is mostly of stellar dynamical nature, our 
approach is grounded on a new computer code for cluster 
dynamics. Its "core" version treats the relaxational evolu- 
tion of spherical clusters. This is the object of the present 
paper. The inclusion of further physical ingredients such 
as stellar collisions and tida l disruptions will be ex plained 
in a complementary ar ticle ( Freitag fc Benz 2001d , see also 
Freitag fc Benz 200 lb , P). In order to obtain a realistic pre- 
scription for the outcome of stellar collisions, we compiled 
a huge database of results from collision simulations per- 
forme d with a Sm oothed Particle Hydrodynamics (SPH) 
code ( |Bcnz 1990|). This work will be described in a third 
paper ( Freitag fc Benz 2001a ). 



This paper is organized as follows. In the next sec- 
tion, we briefly review the available numerical schemes 
that have been applied in the past to simulate the evolu- 
tion of a star cluster and explain the reasons that led us 
to choose one of these methods. Sections ^ to |^ contain 
a detailed description of the code. Test calculations are 
presented in section |^. Finally, in section ^, we summarize 
our work and discuss future improvements. 



2. Choice of a simulation method 

In the previous section, we briefly reviewed the main pro- 
cesses that should play an important role in the evolution 
of galactic nuclei. To treat them realistically, any stellar 
dynamics code has to meet the following specifications. It 
must: 

— Allow for non-isotropic velocity distributions. This is 
especially important to realistically describe the way 
tidal disruptions deplete very elongated orbits (the 
loss-cone phenomenon). 

— Incorporate stars with different masses. Otherwise, 
mass segregation would be neglected and considering 
the outcome of stellar collisions properly would be im- 
possible. Furthermore, it must be possible to make the 
stellar properties (mass, radius) change with time to 
reflect stellar evolution. 

— Be able to evolve the cluster over several relax- 
ation time-scales. The amount of required computing 
("CPU") time can be kept to a reasonable level only 
if individual time steps are a fraction of the relaxation 
time scale rather than the orbital period. 

Several techniques suitable for simulations of cluster 
evolution have been proposed in the literature. Many text- 
books and papers review these methods so we can restrict 



Att anassoula 2000) 



ourselves to a short overview ( 


Binney fc Tremaine 1987 


Spitzer 1987; 


Meylan fc Heggie 1997; Heggie et al. 1999 


Spurzem 1999 


; Spurzem fc Kugel 2000|). 
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As they directly integrate the particles' equations of 
motion, iV-body simulations are conceptually straight- 
forward and do not rely on any important simplifying 
assumptions. Unfortunately, to correctly compute relax- 
ation effects, forces have to be evaluated by direct summa- 
tion (see, however, Jernigan & Porter 1989; McMillan & 
Aarset ~1993| , for possible alternatives). Hence, even with 



specialized hardware such as GRAPE boards (Makino 



1996 |; [iTaiji et al. 1996| ; [Spurzem fc Aarseth 1996| ), fol- 
lowing the trajectory of several 10'^ stars over tens of re- 



laxation times remains impossible (but see Makino 2000) 
for a peek at the future of such systems). A solution 
would be to extrapolate the results of an iV-body sim- 
ulation with a limited number of particles to a higher TV. 



(McMillan 1993 




Heggie 1997 


; Aarseth & Heggie 1998; 


Heggie et al. 199^ 


). The difhculty resides in the fact that 



various processes (relaxation, evaporation, stellar evolu- 
tion . . . ) have time scales and relative importances that 
depend in different ways on the number of simulated stars. 
Thus, in principle, simulating a cluster of stars with a 
lower number of particles, ^ TV*, could only be done 
if all these processes are somehow scaled to the correct 
N^,. Unfortunately, such scalings are difficult to apply, 
precisely because the A^-body method treats gravitation 
in a direct, "microscopic" way with very little room for 
arbitrary adjustments. Furthermore, these modifications 
would rely on the same kind of theory, and hence, of sim- 
plifying assumptions, as other simulation methods. 

To circumvent these difficulties, a class of methods 
has been developed in which a stellar cluster containing a 



very larg e number of bodies is regarded as a fluid ( Larson 
1970b[E; |Lynden-BeU fc Egglcton 1980|; [Louis fc Spurzem 



1991 |; |(piersz fc Spurzem 1994| , |200C| , amongst many oth- 
ers). Such simulations have proved to be very successful 
in disco vering new phenomena like g ravo-thcrmal oscil- 
lations ( Bettwieser fc Sugimoto 1984 ) but rely on many 
strong simplifying assumptions. Most of them are shared 
by Fokker-Planck and Monte Carlo methods (see next 
paragraphs). However, this approach, which is based on 
velocity-moments of the collisional Boltzmann equation, 
further assumes, as a closure relation, some local pre- 
scription for heat conduction through the stellar cluster. 
This is appropriate for a standard gas, where the colli- 
sion mean free path is much smaller than the system's 
size (Hensler et al. 1995). Quite surprisingly, it still seems 



elements (Takahashi 1993, 1995). Here again, the main 
difficulty resides in the treatment of collisions. From a 
practical point of view, these methods represent the clus- 
ter as a small set of continuous distribution functions for 
discrete values of the stellar mass. A realistic modeling of 
collisional effects would then require one to multiply the 
number of these mass classes, at the price of an important 
increase in code complexity and computation time. From 
a theoretical point of view, collisions don't comply with 
the requirement of small orbital changes that is needed to 
derive the FPE. 

So, we have finally turned our attention to the so-called 
"Monte Carlo" (MC) schemes. Even though their under- 
lying hypotheses are similar to those leading to the FPE, 
being particle methods, they inherit some of the A^-body 
philosophy which allows us to extend their realm beyond 
the set of problems to which direct FPE resolutions ap- 
ply. In the Monte Carlo approach, the evolution of the 
(spherically symmetric) cluster is computed by following a 
sample of test-particles that represent spherical star shells. 
They move according to the smooth overall cluster (+BH) 
gravitational potential and are given small orbital pertur- 
bations in order to simulate relaxation. These encounters 
are randomly generated with a probability distribution 
chosen in such a way that they comply with the diffusion 
coefficients appearing in the FPE. 

The most recently invented MC code is the "Cornell" 
program by Shapiro and collaborators (Shapiro 1985, and 
references therein) with which these authors conducted 
seminal studies of the evolution a star cluster hosting 
a central BH. At the end of each time step, the distri- 
bution function (DF) is identified with the distribution 
of test-particles in phase-space. Then the potential is re- 
computed and so are diffusion coefficients that are tabu- 
lated over phase-space. This allows us to evolve the DF 
one step further by applying a new series of perturbations 
to the test-particles' orbits. This code features ingenious 
improvements, such as the ability to follow the orbital mo- 
tion of test-particles threatened by tidal disruption and to 
"clone" test-particles in the central regions in order to im- 
prove the numerical resolution. Despite its demonstrated 
power, we did not adopt this technique because of its al- 
ready important complexity, which would increase signif- 
icantly if additional physics such as stellar collisions, a 
stellar mass spectrum, etc, are introduced. 



to be \alid for a stellar cluster even though the radial ex- 



cursion of a typical star is not 
& Sug: 



■'microscopic" (Bettwieser 



moto 1985; Giersz & Spurzem 1994; Spurzem & 



Spitzer and collaborators (Bpitzer 1975; Spitzer & 
ShuU 19754|b| ; |Spitzcr fc Mathieu 1980D and Henon ( |1966| 



1971b; 1971a; 1973| ) developed simpler MC schemes which 



Takah 4shi 1995 ). However, we discarded such an approach 
because, due to its continuous nature, integrating the ef- 
fects of collisions in it seems difficult. 

A very popular intermediate approach is the "direct" 
numerical resolution of the Fokker-Planck e q uation (FPE) 



(iBinney fc Tremainc 1987| ; [Spitzer 1987| ; [Saslaw 1985| 
for example), either by a finite difference scheme ( Cohn 



1979|,|l|980t |Lec 1987| ; [Murphy et al. 199l| ; [Takahashi 1995| ; 
Drukier et al. 1999| , and many other works) or using finite 



show more potential for our purpose. The simulation of 
relaxation proceeds through 2-body gravitational encoun- 
ters between neighboring shell-like particles. The deflec- 
tions are tailored to lead to the same diffusion of orbits 
that a great number of very small interactions would cause 
in a real cluster. After the encounter, the shell is placed 
at some radius R on its modified orbit. At that time, 
the potential is updated so it remains consistent with the 
positions of the shells. The main difference between the 
"Princeton" code by Spitzer et al. and Hcnon's algorithm 
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is that the former uses a fraction of the orbital time as 
its time step while the latter uses a fraction of the relax- 
ation time. It follows that out-of-cquilibrium dynamical 
processes, like violent relaxation, can be computed with 
the Princeton code whereas Henon's code can only be ap- 
plied to systems in virial equilibrium, which imposes an 
age greatly in excess of the relaxation time. Needless to 
say, this restriction is rewarded by an important gain in 
speed. In this scheme, instead of being determined by the 
equations of orbital motion, the radius i? of a given shell 
is randomly chosen according to a probability distribution 
that measures the time spent at each radius on its orbit: 
dP/dR oc l/urad(^), where Vrad{R) is its radial velocity at 
R. i?-dcpendcnt time steps can be used to track the huge 
variation of the relaxation time from the cluster's centre 
to its outskirts. 

We chose to follow the approach developed by Henon 
to write our own Monte Carlo code. It indeed appeared as 
an optimal compromise in terms of physical realism and 
computational speed. On the one hand, it allows for all 
the key physical ingredients listed at the beginning of this 
section. On the other hand, high resolution simulations 
are carried out in a few hours to few days on a standard 
personal computer. This will enable a wide exploration of 
the parameter space. 

Since Henon's work, this approach has been exten- 
sively modified and successfully applied to the study of the 
dynamical evolution of globular clusters by Stodolkiewicz 
(|l98|; [1986| ) and Giersz (|l99|; |2000a| ; |000b|). Another 
Henon-like MC code has recently been written by 



Josh i et al. (|Joshi et al. 2000| ; [Rasio 2000| ; [Watters et al 
2000; oshi et al. 2001). As far as we know, however, no 
one ever applied this simulation method to galactic nuclei. 

The Monte Carlo scheme relies on the central assump- 
tion that the stellar cluster is always in dynamical equilib- 
rium. This is the case for well relaxed systems. Sufficient 
observational resolution has only recently been obtained 
to allow an estimate of the relaxation times in the nearest 



galactic nuclei. Lauer et al. (1998 ) report relaxation times 
Trei of about 7x 10" yrs, 3 x 10^ yrs and 3x 10^ yrs at 0.1 pc 
for M31, M32 and M 33 respectively. As for the Milky 
Way, [Alexander (1999| ) deduces T.^i ~ 3 x 10^ yrs at 0.4 pc 
(core radius), a value that does not change significantly at 
smaller radii if a p oc cusp model is assumed with 

the parametrisation of Genzel et al. (199^ ) for the veloc- 
ity dispersion. Genzel et al. (1994) get ~ 4 x 10'' yrs 
for the central value but the meaning of this value is un- 
clear, as the dynamical influence of the BH was neglected 
in its derivation. These few values indicate that, perhaps, 
only a subset of all galactic nuclei are amenable to the 
kind of approach we are developing. However, these very 
dense environments with relaxation times lower than the 
Hubble time are precisely the ones which expectedly lead 
to non- vanishing rates for the disruptive events we are pri- 
marily interested in. Furthermore, note that, contrary to 
the case of globular clusters, stellar evolution of a popula- 
tion of massive stars (e.g. in a starburst) can probably not 
disrupt dynamical equilibrium (through mass-loss effects) 



as its time scale (> a few 10^ yrs) is longer than orbital 
periods in a galactic nucleus (< a few 10^ yrs). 

3. General considerations 

3.1. Principles underlying code design 

In our work, we focus at the long term evolution of star 
clusters, on time scales much exceeding the dynamical 
(crossing) time, Tdyn — y/R^JiGMci) , where M^i is the 
cluster's total mass and Rd a quantity indicating its size 
(for instance the half-mass radius). We thus make no at- 
tempt at describing evolutionary processes that occur on a 
Tdyn time scale, most noticeably phase mixing and violent 
relaxation, which are thought to rule early life phases of 
stellar systems (see, for instan ce, section 4.7 of Binney fc 
Tremaine 1987 or section 5.5 of Meylan & Heggie 1997 and 
references therein). Hence, we assume that the cluster has 
reached a state of dynamical equilibrium. Its subsequent 
evolution, driven by processes with time scales ^ Tdyn 
(2-body relaxation, collisions, tidal disruptions and stellar 
evolution), passes through a sequence of such states. 

This reasonable restriction allowed Henon to devise a 
simulation scheme whose time step is a fraction of the 
relaxation time instead of the dynamical time. Naturally, 
this leads to an enormous gain in computation speed com- 
pared to codes that resolve orbital processes likeJV-body 
programs or the Princeton Monte Carlo code (Spitzer & 



Hart 1971; Spitzer fc Thuan 1972) 



Another strong simplifying assumption the scheme 
heavily relies on is that of spherical symmetry. This makes 
the cluster's structure effectively one-dimensional which 
allows a simple and efficient r epre sentation for the grav- 
itational potential (see Sect. 5.1) and the stars' orbits 
and furthermore leads to a straightforward determina- 
tion of neighboring particle pairs. Of course, such a ge- 
ometry greatly facilitates the computation of any quan- 
tity describing the cluster's state, such as density, velocity 
dispersion and so on. An obvious drawback of this as- 
sumption is that it forbids the proper treatment of any 



non-spherical feature as overall rotation (Arabadjis 1997 



Einsel 1998; Einsel & Spurzem 1999|; Spurzem & Einsel 


1999 


; Boily 2000), an oscillating or a binary black hole 


(Begelman et al. 1980; Lin & Tremaine 1980; 


Makino 1997; 


Quinlan & Hernquist 1997; Gould & Rix 2000 


;|Merrit et al. 


2000 


) or an accretion disk interacting with the star clus- 


ter ( 


Syer ct al. 1991; Ranch 1995; Armitage et al. 1996; 


Vokrouhlicky & Karas 1998a, b). 



In Henon's scheme, the numerical realization of the 
cluster is a set of spherical thin shells, each of which is 
given a mass Af , a radius R, a specific angular momen- 
tum J and a specific kinetic energy T. In the remainder of 
the paper, we refer to these particles as "super-stars" be- 
cause we regard them as spherical layers of synchronized 
stars that share the same stellar properties, orbital pa- 
rameters and orbital phase, whose mass is spread over the 
shell and which experience the same processes (relaxation, 
collision, ... ) at the same time. In recent implementa- 
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tions of Hcnon's method by Giersz (1998) and Joshi ct al 
(2000|)] each particle represents only one star. This avoids 
scaling problems connected with the computation of the 
rate of 2- (or many) body processes but would impose too 
large a number of particles for galactic nuclei simulations 
(10^ — 10^ stars). In our code, each super-star consists of 
many stars. Hence, a cluster containing stars may be 
represented by a smaller number of super-stars, A^ss < -^*- 
The number of stars per super-star, N^/N^s, is the same 
for every super-star, in order to conserve energy and an- 
gular momentum (as well as mass when collisions are in- 
cluded) when simulating 2-body processes. 

The super-stars' radii being known, the potential can 
be computed at any time and any place so that the orbital 
energies of all super-stars are straightforwardly deduced 
from their kinetic energies and positions. Hence the set of 
super-stars can be regarded as a discretized representation 
of the one-particle distribution function (DF) / which, 
by virtue of Jean's theorem, only depends on the specific 
orbital energy E and angular momentum J: 

1 



fix, v) = FiEix, v), Jix, v)) = F($s(i?) 



where x and v are the position and velocity of a star, 
R = \\x\\, V — \\v\\ and Vtg is the tangential component 
of V. However, whereas a functional expression of the DF, 
although a complete description of the stellar systerrj^, 
would impose lengthy integrations (resolution of Poisson 
equation, as needed in direct FP methods) to yield the 
gravitational potential, the Monte Carlo representation of 
the cluster provides it directly. From this point of view, the 
Monte Carlo method is closer to an A^-body philosophy 
than to direct FP methods. 

The main difference between the MC code and a spher- 
ical ID iV-body simulation is that the former does not 
explicitly follow the continuous orbital motion of particles 
which preserves E and J. However these orbital constants, 
as well as other properties of the super-stars, are modi- 
fied by "collisional" processes to be incorporated explic- 
itly, namely 2-body relaxation, stellar collisions and tidal 
disruptions. So, in the version of the code described here, 
the MC simulation proceeds through millions to billions 
of steps, eac h of t hem consisting of the selection of super- 
stars (Sect. 4.2.2 ), the modification of their properties to 
simulate the effects of relaxation (Sect. f4.2[ ) and the selec- 
tion o f ra dial positions corresponding to their new orbits 
(Sect. U). 



3.2. Physical units 

In the rest of this paper, unless otherwise stated, we use 



the 'W -body" units as prescribed by [Heggie fc Mathieu 
1986| ). In this system, 



G — 1 (Gravitational constant), 

Mq — 1 (Total initial cluster mass), 

Eq — —1/4 (Total initial cluster energy). 



^ More precisely, the DF, being continuous, is an idealized 
(or statistical) description of the A^* star system. 



It follows that the corresponding unit of length, U\ and 
unit of time, Ut, can be expressed in any system of units 
as 



^ ^ and - 



-4E, 







3/2 " 



(1) 



If the initial cluster's total gravitational energy is written 
as 



f/n = 



qGMl 
'2 i?h 



(2) 



where q is a dimension-less constant and i?h is the half- 
mass radius, we get, for the time unit. 



,-3/2 I 



GMn 



(3) 



For instance, for the often used Plummer model, with 
stellar density p(i?) oc (l + iR/Rp)^y^^^, the 'W-body" 
units are: 



= ( ^) Rp and Ut - 



3/2 



I Rl 
GMo ■ 



As we study systems that are stationary on orbital time 
scales and whose long-term evolution is driven by relax- 
ation, it is more adequate to adopt a time unit that scales 
with relaxation time rather than dynamical time: 



Ut = 



Ut = 7.25(7-3/2Th 



(4) 



where is the standard "half-mass" relaxation time 



( ISpitzer 1987| , p. 40). 

See next section for further explanations of the relax- 
ation time and, in particular the value of 7. 



4. Relaxation. 

4.1. Summary of standard relaxation theory. 

The theory of relaxation in stellar systems can be found 
in classical textbooks ( Hcnon 1973| ; Saslaw 1985; 3pitzei 



1987; Binney & Tremaine 1987, for instance) and will not 
be presented here. However, the treatment of relaxation 
is the backbone of Henon's Monte Carlo scheme. Hence 
a short summary of these issues is particularly worth- 
while, not only to expose the inner workings of the MC 
method, but also to understand the limitations it suffers 
from (as other statistical cluster dynamics approaches) 
that stem from relaxation theory's simplifying assump- 
tions. Furthermore, its specific advantages are also to be 
explained in that framework. 

The basic idea behind the concept of relaxation is that 
the gravitational potential of a stellar system containing 
a large number of bodies can be described as the sum 
of a dominating smooth contribution ($s) plus a small 
"granular" part ((5$). When only the former is taken into 
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account, the phase-space DF of the cluster obeys the col- 
lisionless Boltzmann equation. In the long run, however, 
the fluctuating makes E and J slowly change and the 
DF evolve. The basic simplifying assumption underlying 
relaxation theory is to treat the effects of (5$ as the sum 
of multiple uncorrelated 2-body hyperbolic gravitational 
encounters with small deviation angles. Under these as- 
sumptions, if a test star "1" travels through a homoge- 
neous field of stars "2" which all share the same proper- 
ties (masses and velocities) during St, its trajectory will 
deviate from the initial direction by an angle 9 with the 
following statistical properties: 



fst 



cluster is grossly inhomogeneous, with large density gra- 
dients. Applying Eq. ^ for realistic clusters forces us into 
assuming the "local approximation" , i.e. stating that typ- 
ically b ^ Rc\. Then, not only can we neglect the effect 
of $s during an encounter and treat the trajectories as 
Keplerian hyperbolae, but, as an added benefit, we can 
use the local properties of the cluster (density and veloc- 
ity distribution) as representative of field stars met by the 
test-particle. Admittedly, this is a bold assumption only 
partially justified by the "ln(62/^i)" argument. The valid- 
ity of these approximations has been assessed by compar- 
ing results of codes based on relaxation theory with N- 



= 0, 


(Giersz & Heggie 1994a; ^purzem I 


si Aarseth 199fj; Giersz 




, , 1998 




Portegies Zwart et al. 1998 




Takahashi & Portegies 


'J:i 87l7ihl — ^ = — OL. 
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Zwart 1998; 3purzem 1999). 





In this equation, n is the number density of stars, M\ 
the mass of the test star, Mi the mass of each field-star, 
t^rei the relative velocity between the test star and the 
field stars, 5o is the impact parameter leading to a devia- 
tion angle of 7r/2 and 6niax is a cut-off parameter needed 
to avoid a logarithmic divergence. This ill-defined value 
represents the largest possible impact parameter and is 
thus expected to be of the order of the size of the stellar 
system, i?ci. If is the velocity dispersion in the clus- 
ter and Af* the average stellar mass, the argument of the 
so-called "Coulomb logarithm" can be approximated by 



G (Ml -I- Ma) GAh 



^^^1 Al 



(6) 



where 7 is a non-dimensional proportionality constant. 
This proportionality only holds for a self-gravitating, viri- 
alized cluster. The exact value of 6max/&o is of little im- 
portance as it is embedded into a logarithm. Hence, in 
most applications, a constant 7 is used whose value is 



determined either from theoretical arguments ( spitzer & 



Hart l! 
& Heg 


)71; 


Hcnon 1975) 


^ie 1994a 




1996; 



al. 19991) . The latter 



approach supports the classical "weak encounters" relax- 
ation theory described here by showing good agreement 
with it for properly fitted 7 values. Furthermore it con- 
firms Henon (1975| ) who derived 7 ~ 0.10 — 0.17 for sin- 
gle mass clusters and demonstrated the need for a much 
smaller value when an extended stellar mass spectrum is 
treated. Here, we use 7 = 0.14 and 7 = 0.01 respectively. 

As Eq. ^ is of central importance for the simulation 
of relaxation in the Monte Carlo code, we comment on 
the main assumptions on which its derivation relies. First, 
as already stated, deflections are assumed to be of very 
small amplitude and uncorrelated with each other. The 
contribution of encounters with impact parameters be- 
tween b\ and 62 is of order In (62/^1) so equal logarith- 
mic intervals of b contribute equally to (0^) and most of 
the relaxation is indeed created by "distant encounters" : 
6 > 60 ^ 6* < 7r/2. 

Moreover, the derivation applies in principle only to 
homogeneous systems with a flnite size. However, a real 



Recasting Eqs. |5|, |6| into 



St 



6t 



T. 



(1,2)^ 



we define 



^(1.2) 



rel 



32 ln(7iV,)G2n(Mi -hMa) 



2 ■ 



(7) 



(8) 



We call this quantity the "encounter relaxation time" to 
insist on its depending on the properties of one particular 
class of encounters, namely those between stars of mass 
Ml and stars of mass M2 of (local) density n with rela- 
tive velocity iij-ci- It can be loosely interpreted as the time 
needed for encounters with stars "2" to gradually defiect 
the direction of motion of star "1" by a RMS angle 7r/2. 

4.2. Monte Carlo simulation of relaxation. 

4.2.1. Elementary numerical encounter. 

Contrary to Fokker-Planck codes, Henon's method was de- 
vised to avoid the computational burden and the necessary 
simplifications connected with the numerical evaluation of 
diffusion coefficients (DCs). It does so through a direct 
use of Eq. ^ whose repeated application to a particular 
super-star "1" is equivalent to a Monte Carlo integration 
of the DCs0, provided the properties of field particles "2" 
are correctly sampled. Under the usual assumption that 
encounters are local, this latter constraint is obeyed if we 
take these properties to be those of the closest neighboring 
super-star. Furthermore, this allows us to actually modify 
the velocities of both super-stars at a time, each one acting 
as a representative from the "field" for the other. Hence, in 
the Henon code (as well as in ours), super-stars are evolved 
in symmetrical pairs. This does not only speed up the sim- 
ulations by a factor ~ 2, but also ensures proper local 
conservation of energy, a feature which turned out to be 
a prerequisite for correct cluster evolution. Unfortunately 

^ As the super-star is moved around on its orbit between 
two numerical encounters, the procedure is best described as 
an implicit evaluation of the orbit-averaged DCs. 
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this pairwise approach also impose heavy constraints on 
the code's structure and (perhaps) abilities as we shall 
show later on. 

So the elementary ingredients in the heart of Hcnon's 
scheme are simulated 2-body gravitational deflections be- 
tween neighboring super-stars. However, instead of being 
direct one-to-one counterparts to real individual encoun- 
ters - which would lead to much too slow a code with a 
(huge) number of computational steps scaling as N^^^.^ - 
these are actually "super-encounters" , devised to statis- 
tically reproduce the cumulative effects of the numerous 
physical deflections taking place in the real system over 
a time span St. Thus, such a numerical encounter has a 
double nature. First, it is computed as a (virtual) 2-body 
gravitational interaction with deflection angle ^se in the 
pair CM frame. But being in charge of representing all 
the (small-angle) deflections that test-star "1" experiences 
during St when meeting field-stars "2" , it also has to obey 
Eqs. 1^, ^. Consequently, 9se has to equate the root mean 
squared cumulative deflection 




(9) 



This double nature of the encounter is reflected in the 
whole MC scheme that can be regarded either, quite ab- 
stractly, as a stochastic algorithm to solve the Fokker- 
Planck equation or, more simply, as some kind of ran- 
domized N-body scheme. The second viewpoint, though 
it might be misleading on certain occasions, is the one we 
usually adopt as it allows more intuitive reasoning. 

We now describe the computation of a particular nu- 
merical encounter. It decomposes into the following steps: 

0. A pair of adjacent super-stars and a time step St ar e 
chosen by a procedure to be explained in Sect. 4.2.2 . 

^( 1 2) 

1. The local density n entering the determination of Tj,^j' 
in Eq. ^ is estimated. 

2. The super-stars' velocities, Vi and V2 are randomly ori- 
ented while respecting the angular momenta Ji — \\Ji\\ 
and specific kinetic energy Ti = ^Vi^ of both super- 
stars. This sets the CM- and relative velocities vcm 
and Droi- The former defines the encounter CM frame 
while the latter allows 6se to be determined through 
Eqs. ^ and ||. 

3. In the CM frame, the orientation of the orbital plane 
is randomly chosen around the direction of v^ci- ^se 
being known, computing the post-encounter velocities 
in the CM frame is trivial. 

4. These velocities are transformed back to the cluster 
frame where they define new J^s and TjS for both 
super-stars. 

To compute the local density of star n{R) required in 
step 1, we build and maintain a radial "Lagrangian" mesh 
each of whose cells typically contains a few tens of super- 
stars.0 The cells' radial limits are known, as well as the 

^ The mesh is entirely rebuilt when the number of super-stars 
in any of its cells deviates from some prescribed interval, typ- 



numbcr of super-star they contain. Hence, an estimate of 
the local number density is easily computed by dividing 
the total number of stars in the cell where the encounter 
takes place by its volume. Frequent updating (after each 
super-star orbital movement) and occasional rebuilding of 
the mesh introduces only a very slight computational over- 
head, most CPU time being spent in binary tree traversals 
during potential and rank computations (see Sect. 5.1). 

The computations in steps 2-4 are described in ap- 
pendix The only physical content of all this procedure 
resides in the determination of 9se- Everything else is a 
matter of elementary frame transformations and correct 
random sampling of free parameters in the Monte Carlo 
spirit. 

4.2.2. Choice of the relaxation pair. Determination of 
the time step. 

With no other physical process than relaxation included, 
each individual step in our algorithm comprises three op- 
erations: 

1. Selection of a pair of neighboring super-stars to be 
evolved. 

2. Modification of their orbital properties {Ei and Ji) 
through a super-encounter, as explained above. 

3. For each super-star, selection of a new position on the 
(£'j', Jj')-orbit, i.e. determination of a new radius Ri for 
the spherical shell. This comprises an adequate update 
of the cluster's potential. 

At that point, the code cycles back to 1, i.e., another pair 
is chosen and another step begins. This crucial selection 
process is presented in this section. 

The choice procedure is mainly constrained by the ne- 
cessity of allowing super-stars to have individual time- 
steps Sti that reflect the enormous variations of the re- 
laxation time between the central and outer parts of the 
stellar cluster. When collisions are included, the dynam- 
ical range of evolution times can be even larger. Unless 
this specification is met, the code's efficiency would be 
very low as the overall St would have to reflect the very 
short central evolution time. One could also be concerned 
by the orbital time exceeding St for a large fraction of 
super-stars, a situation inconsistent with the "orbital av- 
erage" approach implicitly assumed in the Monte Carlo 
scheme. However, this problem is actually nonexistent in 
a purely relaxational system whose evolution - under the 
assumptions made in the standard relaxation theory - is 
independent of Trciax/Tdyn, provided is large enough 
(> 100). 

The other important constraint is the need to evolve 
both super-stars in an interacting pair. If the same time 
step is not used for both super-stars, energy is not con- 
served and a very poor cluster simulation ensues. But ad- 
jacent super-stars only form a pair during a unique inter- 



ically [12, 36]. In that way, the mesh structure always matches 
quite closely the super-star distribution. 
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action and then break apart as each one is attributed a 
new radius. So, momentary neighboring super-stars have 
to be given similar 6ti. This strongly suggests use of local 
time steps, i.e. 6t should be a function of R alone. 

Naturally, the time steps have to be sufficiently smaller 
than the time scale for the physical processes driving the 
cluster evolution, namely the relaxation in the present 
case. Hence, we impose: 

StiR) < StoptiR) = fstt.iiR) (fO) 

where Tj-d is some kind of locally averaged relaxation time 
defined as (see Eq. H): 



Trcl CX 



(11) 



In (7iV,)G2n(M,)2 

and fst = 0.005 - 0.05 typically. 

In Eq. n (star number density), (v^) (average 
squared velocity) and (M*) (average stellar mass) are R- 
dependent properties of the cluster. As the only role of 
Trei is to provide short enough 6tiS, an approximate eval- 
uation of these quantities (using a coarse mesh or a slid- 
ing average) is sufficient. On the other hand, too short 
6tiS would fruitlessly slow down the code and should be 
avoided by considering Stopt in Eq. |l^ not only as an up- 
per limit for the time step but also as an optimal value to 
be approached as closely as possible. 

As the members of a pair arrived at their present po- 
sition at different times but have to leave it at the same 
time, once the super-encounter is performed, imposing the 
same 6t to both super-stars is impossible. So, building on 
the statistical nature of the scheme, instead of trying to 
maintain a super-star at radius R during exactly dt{R), 
we only require the mean waiting time for super-stars at 
R to be St{R). As explained by Hcnon (1973| ), this con- 
straint is fulfilled if the probability for a pair lying at R 
to be selected and evolved (and thus, taken away from 
R) during a time span dt is Pscicc(^) — dt/5t{R). This is 
realised in the following way: 

— As it would be difficult to define and use a selection 
probability Pscicc which is a function of the continuous 
variable R, we define it to depend on the rank i of the 
pair (rank 1 designates the two super-stars that are 
closest to the centre, rank 2 the second and third super- 
stars, in ascending order of R and so on). For a given 
cluster's state, local relaxation times Tioi are computed 
at the radial position of every pair. Rank-depending 
time steps are defined that obey inequality 

sm < fstt,i{R{i)). 

— Normalized selection probabilities are computed, 

Tt ^ f""^' 1 ' ' 

clcc \ '^) 



(12) 



with St 



from which we derive a cumulative probability, 



Qsclcc(^) — ^ ^ -Psclcc (j ) • 



(13) 



(14) 



At each evolution step another super-star pair is ran- 
domly chosen according to Pscioc- To do this, a random 
number Xrand is first generated with uniform probabil- 
ity between and 1 . The pair rank is then determined 
by inversion of Qs 



clcc ■ 



* ~ Qsclcc("''^iand)- 



(15) 



The binary tree (See Sect. pj| ) is traversed twice to 
find the id-numbers of the member super-stars whose 
(momentary) ranks are i and i + 1. 

— The pair is evolved through a super-encounter, as ex- 
plained in Sect. 4.2.1 , for a time step St{i). 

— After a large number of elementary steps, typically 
Q.5iVss, the St{i) and Pseiec(*) are re-computed from 
scratch to reflect the slight modification of the overall 
cluster structure. 

As evaporation, collisions and tidal disruptions re- 
move stars from the cluster, the number of super-stars 
Nss generally decreases between two successive computa- 
tions of the selection probabilities. To avoid this problem, 
dt, Pscicc and Qscioc are actually defined as functions of 
a;rank = i / Nss g]0, 1]. 

For the sake of efficiency, we should manage to get for 
'3scicc ^ function that is quickly evaluated. We explored 
two solutions. We first mimicked Henon's recipe, using 
the functional forms: 



-fsclcc (0 
Qsclcc (^) 

Q-} (x) 

^sclccV / 

St{i) 



{l + C)CNi 



ss 



{i + CNss){t + CNss-l) 
(1 + C)i 



■CNss 
, CNss 



1 + C 



6t 



-fsclcc (0 



(16) 
(17) 
(18) 
(19) 



The parameter C > controls the probability contrast 
between the centre (short Ti-ei and St) the outer regions 
(long Tj-ci and 6t). The lower the value of C, the more 
centrally peaked is Psoioc. We determine C by least-square 

fitting (3solec(i) to Qrelit) OC J2]=l ^rc/(j)Q- 

However, being rather arbitrary, Henon's parameter- 
ization could lead to values of St that poorly match the 
shape of T].ei(i?) with the effect of forcing low St and hence 
slowing down the simulation. This concern motivated an- 
other approach. The full rank range is sliced in a few (typ- 
ically 20) bins, in each of which the selection probability 
is set to a constant (either proportional to the maximum 
of in the bin or to its mean value). Such a piecewise 
approximation is naturally expected to adjust better to 
any cluster structure, as shown in Fig. |l]. 

An additional constraint to be mentioned is the need 
to restrict the ratio of the longest time step to the shortest 

* In the least square procedure, points are weighted by 
T~^{j) in order to get a better agreement for the most fre- 
quently selected ranks 
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one, (5tmax/<5imin, to ensure that the outermost shells (the 
cluster's halo) evolves correctly. Otherwise these super- 
stars, most probably placed near their apo-centre posi- 
tions, where relaxation (and, hence selection probability) 
is vanishing, would never be given any opportunity to 
visit more central regions where they can react to the 
adiabatic relaxation-induced modification of the central 
potential and experience 2-body encounters^. Finally, for 
reasons to be presented in appendix |b|, it is necessary that 
the selection probability is a decreasing function of R (i.e. 
of the rank). In practice, these added constraints are ap- 
plied to modify unnormalised selection probabilities which 
are then rescaled to 1. Once the probabilities have been 
worked out, we ensure inequality |l^ is satisfied everywhere 
by setting 

St = fst ma.x{Trc\(i) -^sclcc (20) 

where the maximum is taken over all super-stars. 

As super-stars are randomly chosen to be evolved and 
advanced in time, strict synchronization is never realized 
(except at t = 0) . Each super-star k has its own individual 
time t'-'^-* and synchronization is only achieved statistically 
by requiring that, at every stage in the cluster's evolution, 
the expectation values E{t'^^^) of all t'^^'^s are the same. An 
equivalent statement is to impose an equal expectation 
value -Estop ((^i^*^-*) for the individual time increase of any 

^ If super-stars could be attributed individual time steps, 
setting 5t~^ to an orbital average of T~J^ would naturally pre- 
vent super-stars with small enough i?p from "freezing" in the 
cluster's outskirts. Our procedure amounts to such an orbital 
averaging, in a Monte Carlo fashion. It will fail if the time lag 
between two successive selections of a given super-star is not 
small enough compared with the time over which substantial 
alterations of the cluster's structure occur. 




Li . . . . 1 . . . — . 1 . J 0.09 I- J 

10' 2x10' 0.01 0.1 1 

Step Number Time 

Fig. 2. Evolution of the time dispersion of super-stars 
during the simulation of a core-collapsing cluster. Left 
panel: Time versus step number. The median time tmed 
is shown as a solid line. Dotted lines depict the interval 
[imod — tilled + At_|_] comprising 2/3 of super-stars. 
Right Panel: time evolution of the relative time "disper- 
sion" At/ti„od with M = 0.5(At_ + At+). 



super-star during any evolution step. At the beginning of 
a given step, the super-stars are ranked according to their 
distance to the centre. Selection probability and time step 
depend only on the rank number i{k) so that 

Estep (<5t('^^) = ^^clccW X St(^ (21) 

Probability for Time step for 

selecting the this rank, 

super-star with 
rank i. 
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and Eq. ^ yields: 



step 



sclcc 



(0 



5t 



sclc 



= St Vfc 



(22) 



as desired. Figure ^illustrates how the dispersion of super- 
stars' times evolves during a typical simulation. We define 
the global, "cluster" time to be the median value of the 
super-stars' times. 

5. Orbital displacements and potential updating. 

5.1. Potential representation 



In Sect. 4.1, we explain how relaxation theory, as adopted 
in this work, relies on the assumption that the cluster's 
gravitational potential $ can be described as a dominat- 
ing smooth contribution <i>s whose evolution time scale is 
much longer than the typical orbital time plus a relatively 
small fluctuating 5$. This latter contribution being fur- 
ther reduced to the sum of numerous 2-body encounters, 
we are left with the numerical representation of $s. 

As spherical symmetry introduces many simplifica- 
tions, going beyond this central approximation deeply 
built into Henon's scheme, seems nearly impossible. Its 
most prominent merit is to ensure that stellar orbits, when 
considered on time scales much shorter that Jlciax, are 
easily dealt-with planar rosettes. Therefore, angular 
fluctuations are removed by construction and we represent 
$s as the sum of the contributions of the super-stars, i.e., 
spherical infinitely thin shells of stars. As a consequence, 
radial graininess is still present but its effect turns out 
to be insignificant compared to "genuine" relaxation^. To 
support this claim, we switched off simulated physical re- 
laxation and got a cluster showing no sign of evolution 
(apart from Monte Carlo noise) for a number of steps at 
least three times larger than needed to accomplish deep 
core collapse when relaxation is included. 

Between two successive super-stars of rank i and i + 1 , 
the smooth potential felt by a thin shell of mass M with 
radius R G [Ri, Ri+i] is then simply 



R 



(23) 



with Ai = Mbh + ^ Mj and Bi^^ — 



R. 



where Mj and Rj are the mass and radius of the super- 
star of rank j and A/bh is the mass of the central BH. The 
term 0.5M/R is due to shell self-gravitation. To lighten 
notations, from this point on, $s will simply be referred 
to as the "potential" and the symbol $ is re-attributed to 
it. 

At each step in the numerical simulation, two super- 
stars are evolved which are given new radii and masses (if 

® Moreover, unlike physical relaxation which only depends 
on the simulated number of stars A'^, , radial "numerical" re- 
laxation vanishes as we increase resolution (i.e. the number of 
super-stars A^'ss), see appendix^. 



collisions or stellar evolution is simulated). An easy way 
of obtaining exact overall energy conservation and proper 
account of the adiabatic energy drift of super-stars is to 
update the Ai and Bi coefficients after every such orbital 
displacement. Doing so also ensures that we never put a 
super-star at a radius which turns out to be forbidden (ei- 
ther lower than peri-centre or larger than apo-centre) in 
the updated potential. To sum it up, this choice spares 
us much trouble connected with a potential that lags be- 



hind the super-star distribution. Stodolkiewicz (1982) and 
piersz (199S| ) describe these difficulties, as well as recipes 
to overcome them in their prog rams. Similar prob lems are 
certainly present in the code of Joshi ct al. (2000 ) as they 
recompute the potential only after all the super-stars have 
been assigned new radii. However, performing potential 
updates only after a large number of super-star moves has 
advantages of its own. In particular, in such a code, the 
computing time should scale linearly with the number of 
super-stars (for a complete cluster evolution). This also 
allows them to develop parallelized versions of the Monte 
Carlo scheme (Joshi ct al. 200C). 

If we implement the AiS and BiS as linear arrays, how- 
ever, a large fraction of their Nss elements would have to 
be modified after each step, so the number of numerical 
operations required to evolve the system to a given physi- 
cal time would scale like A^'gg. This steep dependency could 
be avoided by using a binary tree data structure to store 
the potential (and ranking) information (Sedgewick 1988). 
This essential adaptation was alluded to by Henon himself 



(Hcnon 1973) who never published it though. 

At any given time, each super-star is represented by a 
node in the tree. Each node is connected to (at most) two 
other nodes that we shall call his left and right children, 
which are themselves, when present, the "roots" of their 
father's left and right "sub-trees". The rules underlying 
the tree structure are that all the nodes in the left 
"child-tree" of a given node correspond to super- 
stars with lower radii and all the nodes in its right 
"child-tree" to super-stars with higher radii. If we 
define CT^ and TZTk to be the sets of nodes in the left 
and right child-trees of node k, then 

Rk > Rm yrn G CTk ^24-) 
Rk < Rm Vm e TZTk. 
The spherical potential is represented by SAk and SBk 
coefficients attached to nodes. A third value, Sik, allows 
the determination of the radial rank of any super-star. 
These properties, illustrated in Fig. ^, are defined by 
1 -|- number of nodes in CTk , 



Sik = 
SAk = 

SBk = 



Mk- 

Mk 
Rk 



E 



Mm and 



^ Rm' 



(25) 



Then, the super-star with rank i can be found and its po- 
tential coefficients Ai and Bi can be computed by travers- 
ing the tree from root to the node associated with this par- 
ticle. During the same traversal, the SAkS, SBkS and SikS 
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Fig. 3. Diagram of a binary tree with 12 super-stars show- 
ing node properties that encode potential (bAk, SBk) and 
rank {Sik) information. 

of the visited nodes may be modified in order to account 
for the removal of this super-star. Afte r the super-star has 
been given a new radius (see Sect. 5^), it is re- introduced 
into the tree through another traversal. Hence, the poten- 
tial and rank coefficients reflect always exactly the posi- 
tions (and masses) of all the super-stars. The potential at 
any radius, as well as the peri- and apo-centre distances 
of a given super-star, can be computed by specific tree 
traversal routines. The number of operations involved in 
any of these tree traversals is proportional to the number 
of levels i.e., if the tree is kept reasonably well balanced, 
C'(log2(-/Vss))- From time to time, (typically after every 
super-star has been evolved 2 times on average), the bi- 
nary tree is rebuilt from scratch, in order to keep it well 
balanced and to remove empty nodes. More details about 
the implementation of this binary tree can be found in 
appendix |c[ 



5.2. Selection of a new orbital position 

Let's consider a star with known energy E and angular 
momentum J orbiting in a spherically symmetric potential 
$(i?). Its distance to the centre R oscillates between Rp 
(peri-centre radius) and i?a (apo-centre radius) which are 
the two solutions of the equation 



vf,^^2E-2^R) 



i?2 



0. 



(26) 



During a complete orbitFI the time spent in radius interval 
[R,R + dR] is dt oc v^^^{R)dR so that the probability 
density for finding the star at R at any random time (or at 
a given time but without any knowledge about the orbital 
phase) is 



dPo, 



dR PorhV,,d{R) 



(27) 



where Porb = 2/^" dRv^^ is the orbital period. 

Our Monte Carlo scheme avoids explicit computation 
of the orbital motion. It instead achieves correct statistical 
sampling of the orbit of any given super-star by ensuring 
that the expectation value for the fraction of time spent 
at R complies with Eq. Let the sought-for probability 
of placing the super-star at R G [Rp,Ra] be /piac(-R) = 
dPpiac/d-R. According to Eq. if the super-star is placed 
at R, it will stay there for an average time St / Psc\cc{R) ■ 
Then, combining both relations, the average ratio of times 
spent at two different radii i?i , R2 on the orbit is 



istay(^l)\ _ /plac(^l)^sclcc(^2) 



istay(^2)/ /plac(P2)fsolcc(Pl) 
l'rad(-R2' 



(28) 



l'rad(^l) 

As a result, this imposes 

/piac(^) oc 



as required by Eq. 27. (29) 



-^sclcc (-^) 
Wrad(i?) 



(30) 



In appendix ^ we explain how we implement such a 
probability function in an efficient way. 

6. Other ingredients 

6.1. Evaporation and tidal truncation 

Despite its long history, the theoretical understanding of 
the processes leadi ng stars to escape a stella r cluster is 
still not complete (Binney fc Tremaine 1987 Sect. 8.4.1; 



Meylan fc Heggie 1997| Sect. 7.3; [Heggic 200q) . Even with- 
out considering interaction with binaries, the global pic- 
ture seems a bit confusing. Nevertheless, for an isolated 
cluster, the basic mechanism is obvious to grasp, at a "mi- 
croscopic" description level: a star can escape after it has 
experienced a 2-body encounter which resulted in an en- 
ergy gain large enough to unbound it, i.e., to get to posi- 
tive energy. Much of the confusion about the prediction of 
overall escape rates amounts to figuring out whether rare 
large angle scatterings that are neglected by the standard 
relaxation theory could dominate this rate. Indeed it can 
be argued that, in an isolated cluster, stars that are only 
slightly bound and could be kicked away by weak scat- 
terings populate orbits with huge periods and spend most 
time near the apo-centre where encounters are vanishingly 

^ As orbits are generally not closed curves but rosettes, a 
"complete orbit" is defined as the segment of trajectory be- 
tween successive passages at the peri-centre. 
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rare (Henon 196C). According to that picture, the escape of binaries in globular clusters (see, for instance Hut et al 



rate could not be predicted by the "standard" relaxation 
theory, because individual "not-so-small" angle scatter- 
ings would dominate it. 

If this is true, as the MC treatment of 2-body encoun- 
ters relies on the assumption of small relative changes in 
orbital parameters, the method cannot be expected to give 
reliable results for the escape rate from an isolated clus- 
ter (Henon 1971a). Some numerical solution to that prob- 
lem was introduced by Giersz (1998 ). However, when the 
cluster's initial conditions are set to represent a galac- 
tic nucleus, the fraction of stars that evaporate during 
10^° years is very small, so that a precise account of this 
phenomenon is not really required. It would anyway not 
make much sense to devise a complicated treatment of 
evaporation from the nucleus while we neglect the inverse 
process, i.e. the capture of stars from the galactic bulge. 
Our procedure is simply to remove any star whose energy 
is positive after a relaxation/collision process. As can be 
seen in Fig. 0, this simple prescription leads to an amount 



1992, for a review). On the other hand, not much has been 



done concerning galactic nuclei (see Gerhard 1994). 

From a dynamical point of view, only hard binaries, i.e. 
star couples whose orbital velocity Worb is larger than the 
velocity dispersion ay in the cluster, have to be considered. 
In dense systems, they act as a heat source by giving up 
orbital energy and contract (thus hardening further) dur- 
ing interactions with other stars. Of course, the fraction 
of primordial binaries to be labeled as "hard" is a priori 
much higher in globular clusters {a^ of order lOkms^"'^) 
than in galactic nuclei {a^ > lOOkms"^). Binaries can 
also be formed dynamically, either by tidal energy dissi- 
pation during a close 2-body encounter ("tidal binaries") 
or as the result of the gravitational interaction between 3 
stars ( "3-body binaries" ) . The cross section for forming a 
tidal binary strongly decreases with relative velocity (at 
infinity) ( Kim fc Lee 1999 ), so that, in galactic nuclei, 
such processes imply hydrodynamic contact interactions 



that are likely to result in mergers (Lee fc Nelson 1988 



of evaporation in good agreement with the result of Giersz |Benz fc Hills 1992| ; |Lai et al. 1993[ ) . Thus these events 
(199^ 



However, for the sake of comparison with globular clus- 
ter simulations, we also introduced the effects of an exter- 
nal (galactic) tidal field. Due to the sphericity constraint, 
the th ree dimensional nature of such a perturbation can- 



are implicitly treated in our code as a subset of all the 
collisions (Freitag fc Benz 2001d). An interesting counter- 
example to which these arguments do not apply is the 
nucleus of the nearby spiral galaxy M33 whose central ve- 
locity dispersion is as low as ~ 20kms^^ (Lauer et al 



not be [respected. We resort to the usual radial truncation 
approach and consider that a super-star with apo-centre 
radius larger than the so-called tidal radius i?tid imme- 
diately leaves the cluster^. The value of i?tid is about the 
size of the cluster's Roche lobe, _Rtid = c-Rgai(A/ci/-^^gai)^^'^ 
with i?gai being the distance to the centre of the parent 
galaxy whose mass is Afgai and c is of order unity. This 
criterion is clearly a quite unrealistic simplification but we 
do not question it in our work as it is used only for com- 
parison purposes. As the transition from an apo-centre ra- 
dius slightly below i?tid to a value slightly larger does not 
imply large changes in the shape of the orbit, the escape 
rate in this model is certainly dominated by small angle, 
diffusion-like relaxation and must be correctly captured 
by our MC approach. 

6.2. Neglect of binary processes 

The formation, evolution and dynamical role of binaries 
in star clusters are complex and fascinating subjects. An 
impressive number of works have been aimed at the study 



1998) so that tidal binaries should have formed at an ap- 



* Actually, the star would still wander through the cluster 
for a period ~ Porb which could be an appreciable fraction of 
the cluster evolution time scale for low A'",. Neglecting that 
fact could lead to a strong disa greement between A'^-body and 



preciable rate (Hernquist et al. 1991). 

The formation rate of 3-body binaries in galactic nu- 
clei is also strongly quenched as compared to globular 
clusters. Indeed, for a self-gravitating system, the total 
number of binaries formed through this mechanism per re- 
l axation time is only of the order of A^3bb ^ 0.1/ (ln A N^, ) 
( iBinney fc Tremaine 198"7| ; [Goodman fc Hut 1993|) and can 
be completely neglected unless evolution, through mass- 
segregation and core collapse, leads to the formation of a 
dense auto-grav itating core containing only a few tens of 
stars ( Lee 1995 ). Finally, another, somewhat exotic, possi- 
bility is the formation of hard binaries by radiative energy 
losses of gravitational waves during close fly-bys between 
two compact stars ( Lee 1993 , for instance). Note that, if 
present, hard binaries would not only have a dynamical 
role but may also destroy giant stars by colliding with 
them dPavies et al. 19981 ). 

For these reasons, we feel justified not to embark 
on the considerable burden that incorporation of binary 
processes in a Monte Carlo scheme would necessitate. 
However, this has been achieved with a high level of real- 
ism by [Stodolkiewicz (198^ ), Giersz (|l99|; pOOOb| ). Such 
a detailed approach is required to obtain reliable rates for 
binary processes of interest, like super-giant destruction 



Fokk er- Planck based models (Takahashi & Portegies Zwart 
199i). In fact, as the star has to find the "exit door" near 
the Lagrange points before it can effectively escape, it may 
stay in the cluster for many dynamical times even though its 
energy is well above the escape energy. Globular cl usters may 
thus contain a large amount of "p otential escapers" ( Fukushige 



fc Heggie 2000 Baumgardt 200C ) 



by encounters with binaries (Davies et al. 1998), but, if 
needed, the basic dynamical effect of binaries as a heat 
source could be accounted for much more easily using the 
same reci pes that proved suitable in direct Fokker-Planck 
methods ( Lee et al. 1991 , for instance). 

It should be noted that, even in the absence of any ex- 
plicit simulation of binary heating, core collapse is anyway 
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halted and reversed in most Monte Carlo simulations of 
globular cluster evolution! This is due to an effect already 
described by I^6non (1975D and 3todolkicwicz (1982): a 
stiff micro-core consisting of one or a few (< 5) super-stars 
becomes self-gravitating and misleadingly mimics a small 
set of hard binaries by contracting and giving up energy 
to other super-stars. Due to the self regulating nature of 
cluster re-expansion (Goodman 1989), this leads to a post- 
collapse evolution of the overall cluster structure that is 
extremely similar to what binaries produce. Unfortunately 
this does not hold true for the very core whose evolution 
(for instance, whe ther it exper iences gravothcrmal oscil- 
lations or not, see Heggie 1994 ) depends on the nature of 
the heat source. 



7. Code testing 

In this section, we briefly describe the results of a series 
of test simulations we conducted in order to check the 
various aspects of our code and the results it produces. 



7.1. Dynamical equilibrium & spurious relaxation 

The most basic test to be passed is to make sure that 
when relaxation and other physical processes are turned 
off, no evolution occurs in a cluster model whose initial 
conditions obey dynamical equilibrium and radial stabil- 
ity. Beyond stability by itself, the main concern is about 
spurious relaxation introduced by the discrete representa- 
tion of the cluster by a set of super-stars. In other words, 
the supposedly "smooth" potential $ still presents radial 
graininess that could induce some kind of unwanted relax- 
ation (Henon 1971a). It can easily be shown that the time 
scale over which the effects of this spurious relaxation may 
become of importance is of order Tgpur ~ fstNssTrdax- 

This effect has been tested in computations presented 
in appendix Their result is that, provided the num- 
ber of super-stars is larger than a few thousand, there is 
no sign of significant spurious evolution after a number of 
numerical steps larger than what is required in any "stan- 
dard" simulation. No relaxation being simulated, these 
bare-bones steps only consist of orbital displacements as 



10= 



described in Sect. 5.2. Consequently, it appears that these 



radial movements do not introduce appreciable spurious 
relaxation and that the orbital sampling proceeds cor- 
rectly. 

7.2. Core collapse of an isolated single mass cluster 

The next-to-simplest step was to plug in 2-body relaxation 
and to find out whether we could reproduce the well stud- 
ied evolution of an isolated star cluster in which all stars 
have the same mass. We chose a Plummer model because 
it has been extensively used in the literature. Previous re- 
sults for the collapse time are reviewed in Table 0. We 
refer to these many references for a description and expla- 
nation of the physics of core collapse and concentrate on 
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Fig. 4. Code benchmarking. CPU time to deep core col- 
lapse ($(0) = —10) is shown as a function of particle 
(super-star) number for simulations of single-mass isolated 
Plummer clusters. Two sets of simulations were run with 
a different resolution for the radial mesh used to evaluate 
the density and the time steps. Open and black diamonds 
come from runs with 25 and 100 super-stars per cell, re- 
spectively. The factor 2 difference in Tcpu between both 
series probably originates in the fact that, due to Eq. ^ 
the mean time step is sensitive to cells with exceptionally 
low values of Tj-oi (and/or Tcoii if collisions are present). 
Such out-lying values are due to the noise in the grid- 
evaluated Troi and are smoothed out when averages are 
computed with higher number of super-stars. The lines 
are Tcpu — ciA^ss logio(c2^ss) relations computed from 
least square adjustments on points for Nss > 2000. CPU 
times for simulations with Nss < 2000 are probably domi- 
nated by input-output and other system operations rather 
than by the MC algorithm itself. 



some diagrams that describe our simulation of this system 
and allow comparisons with others. 

The results shown here are taken from calculations 
with 512 000 and 2 x 10^ super-stars which took slightly 
more than 100 and 500 CPU hours, respectively, to com- 
plete on a PC with a 400 MHz Pentium II processor. 
Benchmarking of the code is presented in Fig. ^ where 
we plot the CPU time required to attain a value of 
$(0) = — 10 for the central potential as a function of the 
number of super-stars used in the calculation. As most 
computing time (Tcpu) is spend in binary tree traversals 
with a number of operations that scales logarithmically, 
we fitted this data with a relation 



Tcpu = ciiVss logio(c2^ss 



(31) 
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Table 1. Various published values for the core collapse time of an isolated single mass Plummer cluster. Times are 
given in units of half- mass relaxation time T^^^ ( ^pitzer 1987| ). For a Plummer model = 0.093 Z^t (Eq. |). Numbers 
in parenthesis in the last column give the number of particles used. In the "method" column, FP stands for direct 
Fokker-Planck resolution and HMC for Henon-like Monte Carlo schemes. 



Reference 



Henon ( |1973| ; |1975|l 



Numerical method Core collapse time 



Spitzer fc S huU ( |l975aD 
Cohn ( |l979| ) 

Marchant & Shapiro (|l980|) 



Cohn (1980) 



Stodolkiewi cz (|l982| ) 
Takahashi ( |1993D 
Giersz & Heggie ( |l994b| ) 
Takahashi ( |l995| ) 
Spurzem fc Aa rseth (199C) 



Makino (|l996D 
Quinlan ( p96| ) 

Giersz (private communication) 
Lee (private co mmu nication) 
Drukier et a l. (|l999|) 
Joshi et al. ( pOOOD 
This work 



HMC ~ 14.0 - 18.3 (Ik) 

Princeton MC ~ 14.0 - 15.4 (Ik) 

Anisotropic FP 15.9 

CorneU MC 14.7 

Isotropic FP 15.7 

HMC ~ 15.7 (1.2k) 

Isotropic FP 15.6 

iV-body ~ 17.4 (2k) 

Anisotropic FP 17.6 

iV-body 18.2 (2k), 18.0 (10k) 

iV-body 16.9 (8k), 18.3 (16k), 17.7 (32k) 

Isotropic FP 15.4 

HMC ~ 18.3 (4k), ~ 17.5 (64k), 17.4 (100k) 

Isotropic FP 16.1 

Anisotropic FP 17.8 (A^, = 8000), 18.1 (A^, = 50 000) 

HMC 15.2 (100k) 

HMC 17.8 (512k) 17.9 (2000k) 
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Fig. 5. Evolution of Lagrangian radii in a core collapse 
simulation of an isolated single mass Plummer model con- 
sisting of 2 X 10® super-stars. Each curve depicts the radius 
of a sphere that contains the fraction of the mass indicated 
on the label at the right end. These fractions are given 
with respect to the remaining cluster's mass which pro- 
gressiv ily decreases due to evaporati on o f stars. "N-body 
units" ,Ui and Ut are used (see Sect. 3^ 



In Fig. 1^, we present the evolution of Lagrangian radii, 
i.e., radii of spheres that contain a given fraction (0.1% 
to 99 %) of the cluster's mass. In Fig. ^, a subset of these 
radii are used in a comparison with a simulation by Mirek 
Giersz (private communication). We also plot the evolu- 
tion of the decreasing total mass. Clearly, the agreement 
is quite satisfactory. The most obvious difference lies in 
our run leading to a somewhat slower evolution, which 
translates into a core collapse time T^c larger by 2 — 3 %. 
Given the considerable dispersion present in the literature 
for the value of T^c, we judge this discrepancy to be only 
of minor importance. First, due to the stochastic nature 
of Monte Carlo simulations, various runs with the same 
code but using different random sequences yield results 
that differ slightly from each other. This effect is proba- 
bly of minor importance for particle numbers as high as 
2 X 10® but may affect Giersz's data.^ More importantly, 
we stress that although it also stems from Henon's scheme, 
Giersz's code inherited the deep modifications proposed by 
Stodolkiewicz and is actually very different from our im- 
plementation. Also, as Giersz thoroughly and successfully 
tested his program against A^-body data, this comparison 
is very valuable in assessing the quality of our own code. 

Figure |^ is a more direct representation of the same 
information. It shows the density profile p{R) at succes- 
sive evolution phases, deeper a nd deeper in the collapse. 
According to semi- analytical (Lyndcn-Bell & Eggleton 



1980; Heggie & Stevenson 1988, for instance) and numer- 
ical (|Cohn 1980|; [Louis fc Spurzem 1991^ [Pakahashi 1995 



Joshi et al. 2000, amongst others) computations, the cen- 



with constant ci and C2. This is to be contrasted with 
direct A^-body integration which, in its simplest form, re- 
quires 0{N^) operations per relaxation time. 



' To check this, we performed 5 simulations with 10^ parti- 
cles, using different random sequences. We got a dispersion of 
only 1% for the value of Tec. 
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Fig. 7. Evolution of the density profile during core collapse for the same simulation as Figs. || and ^. The dashed line 
shows the slope of the power law corresponding to the self-similar deep collapse solution of the Fokker-Planck equation 
(iHcggic fc Stevenson 1988|). 



tral parts of the cluster evolve self-similarly during the 
late phase of core collapse according to a power-law den- 
sity profile p oc R^^ with ^ ~ 2.2, that extends inwards. 
To illustrate and confirm this behavior, in Fig. H, we plot 
d ln(/9) / d l n(_R) versus ln(_R) and compare with curves ob- 
tained by Takahashi (1995| ) with an anisotropic Fokker- 
Planck code. Although the noisy aspect of our Monte 
Carlo data expectedly contrast with the smooth curves 



from Takahashi's finite-difference code, the agreement is 
clear. The progressive development of a R'^ "^ 
ulation is thus well established. 



m our sim- 



As further evidence for the good performance of our 
code, in Fig. ^ we follow the increase of central den- 
sity p{Q) and central 3D velocity dispersion (w^(0)) in 
our model and compare them with data from an isotropic 



Fokker-Planck computation by H.M. Lee (private com- 
munication) . The collapse time of the two simulations be- 
ing quite different (Tec = 17.9 and 16.1 Wt, respectively), 
we rescaled the time scale by T~^^ in these diagrams to 
get more meaningful comparisons. Here again, the level 
of agreement is more than satisfactory. Incidentally, we 
note that a (w^(0)) oc p{0)'^ relation is quickly estab lished 
with C ~ 0.10 as previously noted by many authors (Cohn 
1979|, |1980|; [Marchant fc Shapiro 1980|; [Takahashi 1995|). 



Finally, in Fig. [T^, we investigate the evolution of 
anisotropy, measured by the usual parameter A — 2 ~ 
where wtg and Wiad are the tangential and 
radial components of the stellar velocity. The develop- 
ment of a radial anisotropy at every Lagrangian radius 
is clearly visible. Our curves are very similar to those 
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Fig. 6. Comparison of our collapse simulation (solid lines) 
with a computation by Mirek Giersz (dash-dotted lines). 
Our simulation is the same as in Fig. |. Giersz used 100000 
super-stars. The top panel shows the Lagrangian radii. 
The bottom panel depicts the total mass. Contrary to 
ours, Giersz's code is able to go past core collapse by sim- 
ulating the formation of 3-body binaries and their giving 
up energy to other stars. 



obtained by Takahashi (1996 ) and Drukier et al. (1999 ) 
using anisotropic FP codes. Globally, although the agree- 
ment is not as close as in previous diagrams, this rela- 
tive mismatch is weakened when examined in the light 
of the differences between both FP simulations. It thus 
seems reasonable to think that these differences could be 
due to the further simplifying assumptions required in the 
derivation of direct anisotropic FP schemes. We also in- 
clude the Monte Carlo simulation by Giersz (199^ ) in the 
comparison. Due to the relatively low number of particles 
used by this author (10^), this data contains large sta- 
tistical fluctuations. It is nonetheless clear that it shows 
significantly more radial anisotropy in the central regions 
at early times, compared to the other simulations. The 
reason for this difference is unknown to us but may be 
due to the use of radial "super-zones" by Giersz to define 
block time-steps. In his scheme, super-stars from an in- 
ner super-zone are not allowed to skip to an outer zone, 
unless both zones happen to be synchronized at the end 
of their respective time-steps. This clearly could lead to 
some "particle restraining" in the central parts but it is 
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Fig. 8. Logarithmic density gradient for successive stages 
of the core collapse of the Plummer model. Points rep- 
resent our data for the same times as in Fig. |^ (except 
the initial state which is not represented here). Curves 
arc from a simulation by Takahashi (1995) for stages in 
core collapse with about the same central density in- 
crease. The leftmost curve corresponds to a collapse phase 
slightly deeper than attained by the Monte Carlo simula- 
tion. Dotted vertical lines show the decreasing values of 
the core radius = ^/3{v^{0)) / {ATrp{0)). 



unclear why this effect would appear in the anisotropy pro- 
files but not in the density data. Comparison with data 
from an A^-body code would allow us to settle these ques- 
tions. Unfortunately, present-day accessible N values (a 
few 10^) yield anisotropy curves still too noisy to be of 
real use ( Spurzem fc Aarseth 1996| , for instance). 

To summarize this subsection, we can safely conclude 
that, when applied to the highly idealized relaxation- 
driven evolution of a single-mass cluster, our Monte Carlo 
implementation produces results in very nice agreement 
with many other modern numerical methods and theoret- 
ical predictions. To step closer to physical realism, we now 
present the results for multi-mass models. 

7.3. Evolution of clusters with two mass components 

Cluster models with stars of identical masses fall short of 
any realistic description of real clusters. Indeed, it has long 
been known that the evolution is profoundly affected by a 



stellar mass spectrum (inagaki fc Wiyanto 1984 ; Inagak 



|fc Saslaw 19851 ; [Murphy fc Cohn 1988| ; [Takahashi 1997| , for 
instance). If the cluster is initiated with the same spatial 
and velocity distributions for the various stellar masses, 
2-body gravitational encounters will attempt to enforce 
equipartition of kinetic energy between stars of different 
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Fig. 9. a) Evolution of the central density during the 
core collapse of the Plummer Model. On this diagram, the 
line for our simulation (2 x 10^ particles) cannot be dis- 
tinguished from an isotropic Fokker-Planck model by Lee 
(private communication)! The time scales of both simu- 
lations have been scaled by their respective core-collapse 
times. A slight smoothing has been applied to our data, 
b) Same as panel a), but for the 3D central velocity dis- 
persion. The solid curve shows our simulation, the dash- 
dotted line is the model by Lee. 



masses, hence causing heavy stars to segregate towards the 
centre. Depending on the relative masses and numbers of 
heavy and light st ars, a temporary equipartition may be 



(nearly) attained ( Inagaki fc Wiyanto 1984 ; Watters et al 



2000|). [The latter case corre sponds to the segregation (or 
"stratification") instability ( Spitzer 1987 ) but even when 
it does not occur, a central sub-cluster containing massive 
stars forms and collapses quickly because its relaxation 
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Fig. 10. Evolution of the average anisotropy parameter 
within shells bracketed by Lagrangian spheres containing 
15 to 20%, 45 to 50%, 70 to 75% and 90 to 95% of the 
cluster mass. Solid lines show the data from our simulation 
with 2 X 10^ particles. Short dash curves are results from 
Takahashi (l"996| ) for the anisotropy at 20%, 50%, 75% and 
90% Lagrangian radii. Long dash curves are taken from 
Drukier et al. (1999D for 50%, 70% and 90% radii. Dot- 
dashed lines are the (smoothed) results of a 10^ particle 



simulation by [Gicrsz (1998D for 14-20%, 44-50%, 70-76% 
and 90-96% shells. The time scales have been scaled by 
their respective core-collapse times. A slight smoothing 
has been applied to our and Giersz's data. 



time is much lower than the overall value. For a two com- 
ponent model, with stellar masses mi and m2 {m2 > mi ) , 
the time scale for equipartition and induced segregation 



can be as low as Toq ~ toi/to2 • Jlci ( Spitzer 1969| ). As a 



consequence, the structure of multi-mass clusters evolves 
much faster than single mass models. 

As a first validation of our code in the multi-mass 
regime, we simulated clusters with two mass components. 
Such models have been used by p^jce (1995 ) to study the 
fate of stellar black holes (SBHs) in galactic nuclei. He as- 
sumed a simple stellar population with all main sequence 
stars with mass mi = 0.7 Mq and a given fraction of 
7712 = 10 Mq SBHs. For a board range in the total mass 
fraction of SBHs, his 1-D Fokker-Planck simulations show 
a collapse time for the SBH sub-system that is reduced by 
factors of tens compared to the single-mass case. As shown 
in Fig. 11, we successfully reproduce his results for various 
mass fraction of SBHs. Note that the MC method is un- 
able to simulate reliably clusters containing a very small 
SBH fraction if the number of super-stars has to be much 
lower than the number of simulated stars. For instance, if 
we wanted to simulate a cluster with MsBn/Mtot — 10""^ 
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Fig. 11. Collapse times for cluster models consisting of 
2 components: MS stars with mass mi and stellar BHs 
with mass TO2. The initial clusters are Plummer mod- 
els without segregation. The collapse time Tec is nor- 
malized by the collapse time for a single-mass cluster, 
ric'™ ''- Collapse times are plotted as a function of the 
mass fraction of SBHs in the cluster. Black dots show 
our MC simulations for mi = 0.7 Mq and = 10 Mq. 
The leftmost point actually comes from a simulation with 
-^SBn/Aitot = 0. The circled dots are for simulations with 
512 000 super-stars; 64 000 super-stars have been used 
for other simulations. The open diamonds connected by 
the dashed line are data from 1-D Fokker-Planck simula- 
tions by Lee (1995| ). The open circles show our results for 



m = 0.3 Mq and M = 10 with 512 000 super-stars. 



where Msbh is the total mass of SBHs, with 5 x 10^ super- 
stars, only 35 of them would represent SBHs, a number 
clearly too small to resolve the collapse of the central sub- 
cluster of SBHs. To illustrate the process of mass segre- 
gation, in Fig. ^ we plot the evolution of the Lagrangian 
radii of both mass components for two of these cluster 
models. 

7.4. Evolution of a tidally truncated multi-mass cluster 

To study the evolution of a cluster with a continuous 
mass function, we simulated a model with initial con- 
ditions set according to Heggie's "collaborative experi- 
ment "[] ( iHcggie et al. 19991). This is a King model with 
Wo = -<l>(o)/g^ = 3 dBinney fc Tremaine 1987| , p.232). A 
spherical tidal truncation is imposed at i?tid — 30 pc. The 




Time in 10^ years 

AfsBH/Mtot = 0.81 




Time in 10^ years 

Fig. 12. Evolution of the Lagrangian radii of the compo- 
nents for two models of Fig. |ll|. The solid and dashed lines 
are for stars with mi = 0.7 Mq ("main sequence": MS) 
and m2 = 10 Mq ("stellar black holes": BH), respectively. 



mass spectrum is: 

diV 



M-2-35 for 0.1 Mq < M* < 1.5 Mp 



All experiment data is available at http:// 
www .maths . ed. ac .uk/~douglas/experiment .html 



and the total mass is 6 x 10'*A/q. Hence, the number of 
stars is N^, = 2.474 x 10^. There is no initial mass segrega- 
tion and no primordial binaries. According to the rules of 
the experiment, no stellar evolution has to be simulated 
but the heating effect of binaries could be incorporated to 
simulate the post-collapse evolution up to complete evap- 
oration of the cluster. 

Our code lacks the ability to simulate the formation 
of binaries and their heating effect. However, as explained 
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Fig. 13. Evolution of the Lagrangian radii for a cluster 
with initial conditions according to Heggie's "collaborative 
experiment" . Solid lines are the result of our simulation 
with 256 000 super-stars. Dashed lines with black dots are 
from a computation by Mirek Giersz (10, 50 and 100% 
radii). The tidal radius is label as "100%" . "TV-body" units 



are used. Unit of time: 
9.55 pc. 



5.70 X lO'^yrs. Unit of length: 



in Sect. 6.2, these processes do not switch on before the 
core has collapsed down to a few tens of stars. As a con- 
sequence, we should be able to tackle the evolution of this 
system up to deep core collapse. 

Many researchers, using a variety of simulation meth- 
ods, from gas models to A^-body codes, have taken part in 
the "collaborative experiment" . Their results show a very 
important dispersion. For instance, the obtained core- 
collapse times range from 9 to more than 14 Gyrs while the 
values for cluster's mass at this time lie between 2.2 x lO'' 
and 4.75 x 10** Mq! In our simulations with 16 000 to 
256000 super-stars, we find a collapse time of 12.5 to 
13.4 Gyrs with a remaining mass varying from 4.63 x 10^ 
to 4.37 X 10" Mq. 

Factor 2 discrepancies can even been found amongst 
simulations using the same scheme, e.g. iV-body codes. 
There is a clear tendency for iV-body to yield values of Tec 
shorter than those produced by other, statistical, meth- 
ods. Another perplexing fact is that the results of A-body 
simulations do not converge to those of statistical methods 
as A^ increases, contrary to naive expectations. The cause 



of thes e unexpec ted results has been traced by Fukushige 



& Heg gie (2000|) to two combining facts. First, for the 
realistic, non-spherical, tidal potential used in those sim- 
ulations, stars with energies above the escape energy can 
stay in the cluster for many dynamical times before they 
actually leave it. Second, the way the models where initi- 




Time 

Fig. 14. Mass segregation diagram for a cluster with ini- 
tial conditions according to Douglas Heggie's "collabora- 
tive experiment" . Each curve show the evolution of the 
average mass of stars in a sphere that contains a given 
fraction of the total cluster's mass, as indicated by labels 
on the right. Solid lines are the result of our simulation 
with 256 000 super-stars. Dashed lines with black dots are 
from a computation by Mirek Giersz (10, 50 and 100% 
radii). The unit of time is 5.70 x lO^^yrs. 

ated led to clusters containing, from the beginning, a large 
fraction of such "potential escapers" , instead of being in 
equilibrium in the tidal field. 

Given this confusing picture, it seems more sensible 
to compare our results to those produced by a similar 
computatio nal approach. Mirek Giersz applied his Monte 
Carlo code ( [Giersz 1998i ^000b| ) to this system. In Fig. ||, 
we show the evolution of Lagrangian radii for his simula- 
tion (up to binary-induced rebound) and ours. Similarly, 
Fig. |lj compares the evolution of the average mass of 
stars. This latter diagram clearly shows how strong mass 
segregation effects are in multi-mass clusters. The rela- 
tively good agreement to be read from these figures sup- 
ports our code's ability to handle star clusters with a mass 
spectrum. 

8. Conclusions 

8.1. Summary 

In this paper, we have presented a new stellar dynam- 
ics code we have recently developed. It can be seen as a 
Monte Carlo resolution of the Fokker-Planck equation for 
a spherical star cluster. Although stemming from a scheme 
invented by Henon in the early 70 's, it was deemed optimal 
for our planned study of the long-term evolution of dense 
galactic nuclei hosting massive black holes. The main ad- 
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vantages of this kind of approach are a high computational 
efficiency (compared to A'^-body codes), on the one hand, 
and the ability to incorporate many physical effects with 
a high level of realism (compared to direct Fokker-Planck 
resolution or to gas methods), on the other hand. These 
features explain the recent revival of Henon's approach in 



the realm of globular cluster dynamics by Giersz (1998) 
and Joshi et al. (2000). To the best of our knowledge. 



however, we are the first to apply it to galactic nuclei (see 
Freitag fc Benz 2001bp|). 



The version of the code presented here only includes 2- 
body relaxation. Spherically symmetric self-gravitation is 



computed exactly. Arbitrary mass spectrum and velocity 
distribiition, isotropic or not, can be handled without m- 
troducing any extra computational burden. The test com- 
putations we carried out allow us to be confident in the 
way our code simulates the evolution of spherical star clus- 
ters over the long term, as driven by relaxation. 

The computational speed of our code is highly sat- 



this list broadly reflects the probable order of inclusion of 
these effects in our simulations. 

— Stellar collisions. 

— Tidal disruptions. 

— Stellar evolution. 

— Capture of stars by the central BH through emission 
of gravitational radiation ( "GR-captures" ) . As these 
events are a very promising source of gravitational 
waves for the future space-borne interferometer LISA, 
reliable predictions for their rate and characteristics 
are highly desirable even though they are unlikely to 
play a dominant ro le in the BH's growth ( Danzmann 



isfying The evolution of a single-mass globular cluster 



with 512 000 super-stars up to core collapse takes about 
5 CPU days on a 400 MHz Pentium-II processor. This 
can be compared with the three months of computa- 
tion required by Makino (1996| ) to integrate a cluster 
with 32 000 super-stars on a GRAPE-4 computer spe- 
cially designed to compute forces in an A^-body algorithm. 
However, the significance of such a comparison is some- 
what blurred by the fact that Makino integrated the sys- 
tem past core collapse and that the hardware in use is 
so different. Nevertheless, the speed superiority of our 
Monte Carlo scheme over A^-body really lies in a CPU 
time scaling as N ■ log(ciV) instead of iV^~^ (Makino re- 
ports Tcpu oc N'^-'^.). Furthermore, Monte Carlo simula- 
tions do not have to resolve orbital time-scales; their time 
step is a fraction of the relaxation time which is of an or- 
der 10^ times larger for a million-star self-gravitating clus- 
ter. This ensures that Monte Carlo simulations will remain 
competitive in the next few years, even after the advent of 
special -purpose A^-body computers with highly increased 



perfo ri iances like the GRAPE-6 system ([Makino 199S , 
200C). Monte Carlo codes like ours are bound to become 
the tool of choice to explore the dynamics of star clusters 
by allowing investigators to run many simulations with a 
variety of initial conditions and physical processes. Run- 
of-the-mill personal computers are sufficient to get quick 
results without sacrificing too much of the physical real- 
ism. 

8.2. Other physical ingredients and future 
developments 

2-body relaxation is only one amongst the many physical 
processes that are thought to contribute to the long term 
evolution of dense galactic nuclei or are of high interest of 
themselves even without a global impact on the cluster. 
Here, we list the most important of them and comment 
on those not already discussed in Sect. |l|. The order in 



1996| ; [Thorne 1998)) 



— Large scale gas dynamics. Gas is released by stars dur- 
ing their normal evolution (winds, SN explosions, . . . 
) or as a consequence of collisions. Recent 2-D h ydro- 
dynamic simulations by [Williams et al. (1999| ) have 
revealed a variety of behaviours that were not cap- 



tured by previous works ( Bailey 1980 ; David et al 
19874 |fc|)- The determination of the fraction of gas ac- 



creted by the central BH and the fraction that escapes 
the galactic nucleus appears to be a difficult but im- 
portant problem. 

Interplay with outer galactic structures. Contrary to 
globular clusters, the nucleus of a non- interacting 
galaxy is not subject to a strong tidal field. However, 
it is not an isolated cluster. It is embedded in a larger 
structure (bar, bulge, elliptical galaxy) whose gravita- 
tional potential is generally not spherically symmetric 
and with which it can exchange stars and gas. 
Interact ions with binary stars. This has been discussed 



in Sect. 3.2 



Other interactions between the central black hole and 
stars. A number of more or less exotic mechanisms 
have been proposed in the literature, most of them 
as alternate mechanisms to feed the central BH with 
stellar fu el. Amongst others, w e mention tidal capture 
of stars ( Novikov et al. 1992), their i nteraction with 



a central accretion disk (Ranch 1995; Armitage et al 



1996, and others), mass transfer to the BH by a close 



orbiting star ( Hameury et al. 1994 ) and the influence 
of the UV/X-ray flux from the accreting BH on the 
structure and evolution of nearby stars (for instance 
X-ray induced stellar winds, see Voit & Shull 1988). 
Cluster rotation. A few recent observations indicate 
that the centre-most regions in a cluster may present 



substantial amounts of rotation (Genzel et al. 2000 
Gebhardt et al. 200(\). 



The original Henon's code was devised to study ideal- 
ized globular clusters whose evolution is solely driven by 
relaxation. Such models are only remotely connected to 
galactic nuclei. Unfortunately, the processes possibly at 
play in galactic nuclei are so numerous and (for many of 
them) uncertain that fully consistent simulations, incorpo- 
rating all the physics, seem to be beyond reach for many 
years still. Such simulations would look misleadingly re- 
alistic but yield little insight into the importance of each 
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individual process and how it interplays with the others. 
To favor such an understanding, we restrict our discussion, 
for the time being, to a few ingredients that are deemed 
particularly important. We want to get a good insight into 
the "workings" of these simplified models before we add 
more complexity and uncertainties by including further 
physics. Some of these ingredients are very likely to play a 
key role in the evolution of the cluster: tidal disruptions, 
stellar evolution, maybe collisions . . . Other processes, like 
GR-captures, may be too rare to have an noticeable influ- 
ence on the overall dynamics and structure of the system 
but have great observational promise as individual events. 



In the following paper of this series (Freitag & Benz 
2001d^ we'll describe how stellar collisions and tidal dis- 
ruptions are treated. The next effects to which we shall 
turn are stellar evolution, includ ed in a simpl ified way in 
the latest version of the code ( Freitag 2000 ), and GR- 
captures, for which enc ouraging results have already been 
obtained ([Freitag 200l|) . 
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Appendix A: Practical computation of a 
super-encounter 

In this appendix, we detail the vectorial operations corre- 
sponding to steps 2, 3, and 4 of an individual encounter 
between two neighbouring super-stars, as described in 
Sect. 



4.2.1 



Step 2. The situation is depicted on top of Fig. |A.l[ A 
local reference frame, at rest in the cluster, is defined with 
axis Oz pointing in the radial direction from the cluster's 
centre and Vi G Oxz. From the specific angular momenta 
Ji, kinetic energies Ti and distances to centre Ri of the 
super-stars, the moduli of the radial and tangential com- 
ponents of the stars' velocities are deduced: 



We set vf 



,tg\2 



and vl 



rad 
^1 1 



(A.l) 

vfcOs{lf2), 



V2 = V2 sin(<p2) and d| = , with a random value for 

the angle 1^2 (s [0, Itt]) and the sign of i>|. This procedure 
accounts for the freedom in the relative orientation of 
and ^2, thus ensuring a correct sampling of the encounters' 
incoming velocities. 

Step 3. We now switch to the encounter CM refer- 



ence frame (see bottom panel of Fig. A.l) which has the 
same axis orientation as the cluster frame but translates 



http : //obswww . unige . ch/~pf emiige/gravitor/ 
gravitor_e . html 




Fig. A.l. The reference frames used in the computation 
of the 2-body encounter. Vi^2 are the velocities of the stars 
in the local cluster frame, 11^1,2, their velocities in the en- 
counter reference frame, Wrci — Vi ~ V2 is the relative ve- 
locity and vqm the velocity of their centre-of-mass (CM) 
in the cluster frame. See text for more details. 

with velocity vqm = Ait>i 4- \2V2 {K = Mi/ {Mi + M2) 
where Mi,2 are the masses of the stars). We use the no- 
tation 1)1.2 for the stars' velocities in the cluster frame 
and Wi^2 when we express them in the encounter frame. 
Prime (') denotes quantities after the encounter. We build 
an orthonormal vector set {e||, e-^, e^} with ~ v-cd/v-cci 
and e^,es -L Vrci- The orientation of the orbital plane 
spanned by {e||,ej^} is set through a randomly chosen 
angle /3 defining ej_ — cos(/3)e-y -I- sin(/3)e5. In this frame, 
the effect of the gravitational encounter is simply to ro- 
tate the initial velocities (at infinity) Wi ~ A2frci and 
-0)2 = — Aifrei by an angle 0se, so: 

w[ = X2V'^^1, w'2 = -Xiv'^^i 
with v[.^^ = Vrc\ (cos(6'sE)e|| + sin(6'sE)e_L) 

Step 4. Finally, the relevant post-encounter proper- 
ties of the super-stars in the cluster frame are given by 
straightforward formulae: 



(A.2) 



T' 



J' 



1 



(A.3) 
(A.4) 



with v'^ — w'. 
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Fig. B.l. Diagram of the random selection of an orbital 
position using a rejection method with a Keplerian com- 
parison function. A constant Psoicc(^) is assumed. 




Appendix B: Random selection of a new orbital 
radius 

Here, we explain how we determine a new radius R for a 
super -sta r after it has experienced a super-encounter (see 
Sect. U). 

Our goal is to generate random R values whose dis- 
tribution density comply with Eq. The main difficulty 
lies in the fact that /piac(^) is not an easily computable 
function. First, Pscicc really depends on the rank rather 
than on R and has g enerally no simple analytical expres- 
sion (see Sect. 4.2.2|) . Furthermore, Urad(-R) is a function 
of $(i?) which is not known analytically either. Obtaining 
the value of the rank or the potential (through local A 
and B coefficients) at any given R implies a tree traver- 
sal. Given these intricacies, we don't even attempt trans- 
forming a uniform-deviate random variable through the 
inverse function of the cumulative probability and turn 
to a rejection method (Press et al. 1992, Sect. 7. 3). The 



trick is to find a comparison function /comp, a more docile 
probability density which can be made a uniform upper 

bound to /plac(i?) (/comp(i?) > a/plac(i?) Vi? G [i?p,i?a] 

for some constant a). We then proceed by generating a 
random number -Rrand according to /comp and another, 
^rand, with Uniform deviation between and 1. If Xjand < 
a/piac(i?rand)//comp(i?rand), -Rrand is kept. Otherwise it is 
rejected and we try again. The accepted -Rrand values fol- 
low the /piac distribution. Of course the closer the com- 
parison function is to /piac(-R), the fewer rejection steps 
and the more efficient is the method (remember that a 
tree traversal is realised for each /piac(-R) evaluation). 

By construction, Pseiec is constrained to decrease with 
rank/radius, so an easy but inefficient upper bound is its 
value at peri-centre. As for Wrad(-R)~^, we first apphed 
Henon's variable change {R = Rp + {Ra — i?p)s^(3 — 2s)) 
to remove the divergences at peri- and apo-centre but 
failed to find a safe upper bound for the resulting proba- 
bility density of s. We instead use the fact that a shifted 



Fig. B.2. Example of the shape of the probability den- 
sity /piac oc /'scioc(-R)/wrad(-R) for placing a super-star 
on its orbit. Henon's s variable has been used instead 
of R to remove the divergences of v~^^ at peri- and 
apo-centre {R ^ Rp + {R^ - Rp)s'^{3 - 2s)). fpUc is 
shown in solid line along with the simple Keplerian up- 
per bound Pscicc{Rp) / v[^^^\r) (dot-dashed fine) and a 
refined Keplerian upper bound (dashes) with coarse ac- 
count for the important radial decrease of Pscioc- Psoioc was 



computed with the second method described in Sect. 4.2.2 



where is evaluated on a coarse rank grid. This is re- 
sponsible for the saw tooth shape of /piac- Obviously, the 
"raw" Keplerian /comp would mostly generate R values 
near the apo-centre very likely to be rejected due to the 
actual very low probability density for such Rs. 



Keplerian potential $Kcp(P) — C\ / R + C2 equal to $ at 
Rp and i?a is everywhere higher (or equal) in between, so 
that 



1 



< 



.rad(i?) - .(f7)(i?) 

\/ RpRa 



R 



2J ^{R - Rp) (i?a - R) 



(B.l) 



Furthermore, the cumulative probability function for this 
Keplerian bound is analytical: 



FKcpiR) « / 

Jn 



■ dr 



Bp \J{r - Rp) (i?a - r 
arctan 



I 



(B.2) 



. , R — i?p , i?a — i?p 
with X = and e = 



i?a — Rp i?a + i?p 

So, to summarize, the selection of a new orbital radius 
proceeds as follows (see Fig. B.l ): 
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1. Generate a random number Xj-and with [0, l]-uniform 
deviate. 

2. Compute i?i.and = i?p + (^a - Rp)F^^^{Xri,nd) using, 
for instance, Newton-Raplison algorithm. 

3. Keep R = i?rand with probabihty 



-fkeep — 



Psolec(i?) v[^r\R) ^ 

-^sclcc 



This procedure actually had to be improved, for Pgeiec 
is generally very rapidly decreasing with rank. As a result, 
if the super-star's orbit spans a large i?-range, a constant 
bound often proves to be highly inefficient, requiring hun- 
dreds of rejection steps. So, when the number of unsuc- 
cessful tries reaches a limiting value, the [-Rp, i?a] interval 
is sliced into a few sub-intervals and a stepped bound on 
Pseiec is devised by computing its value at the lower rank 
of each sub-interval. Hence a comparison function is ob- 
tained that lies closer to /piac- However, in case a large 
number of rejections still fails to select a R value, further 
successive refining is realized, using more and more sub- 
intervals to get closer and closer upper bounds to Pscicc- 
In practice, we use successively 5, 20 and 50 sub-intervals 
when the number of unsuccessful rejection cycles exceeds 
10, 100 and 1000 respectively. This modification clearly 
complicates the computation of -Fkep and its inversion dur- 
ing phase 2. It appears as a numerical overhead as it im- 
poses a tree traversals to determine lower rank values for 
each sub-interval. However, such added intricacies allow 
to break off from a (nearly) never ending rejection cycle. 
The improvement on the Kepleria n com parison function is 
depicted for a typical case in Fig. B^. With this method, 
the average number of rejection cycles to attribute a new 
i? to a super-star lies between 5 and 10 in our cluster 
simulations. 



Appendix C: Binary tree structure and algorithms 

In this appendix, we present in some detail the implemen- 
tation of the binary tree in charge of storing the potential 
and ranking information of the super-star cluster. 

The logical structure of the tree is implemented by 
three integer arrays: l_son, r_son and father. l_son(k) 
is the number of the "left son" node of node k and so on, 



see Fig. C.l. When the tree is (re-)built, each node k is at- 
tributed a super-star super_star (k) . When this particle 
is evolved and moved to another radius, the node becomes 
empty and another leaf node is added to host the super- 
star. Leaving empty nodes simplifies the tree update pro- 
cedures with the cost of increasing its size. This introduces 
some numerical overhead as it causes a faster increase of 
the number of hierarchical levels (particularly in the cen- 
tral regions where the time steps are shorter, see Fig. |C.l[) 
but this inconvenience is probably not a serious concern 
for the tree is rebuilt from scratch from time to time in 
order to keep it reasonably well balanced. This operation 
is computationally cheap compared to the numerous tree 



traversals and would be called for even if empty nodes 
were not created]^ 

To find a super-star with rank i and compute the po- 
tential at its position, one traverses the tree from the root 
to the corresponding node using the algorithm sketched 
in Fig. p.2| a. At the end of this tree traversal, fcscarch 
points to the proper node and the super-star's potential 
is $ = ^^/^fcsoarch — B. A very similar routine is used 
to compute the potential ^{R) at any arbitrary radius R. 
Flow charts for addition/removal of a super-star into/from 
the binary tree are depicted in Fig. C.2b, c. 

Another important role of the binary tree is the com- 
putation of peri- and apo-centre radii for a given super- 
star, an obvious prerequi site to the radial placing pro- 
cedure described in Sect. 5^. Note that a high level of 
precision is called for: unphysical cluster evolution could 
expectedly arise if super-stars' orbits suffered from any 
systematic bias. For lack of explicit knowledge of ^{R), 
solving Eq. |2^ is not straightforward. The main opera- 
tion, depicted in Fig. |C.2|d, is a tree traversal in a search 



for the node kp whose radius lies immediately below the 
peri-centre, i.e. with Q{R) — v'^^^{R) < and dQ/di? > 0. 
This provides us with the local Ap and Bp potential co- 
efficients. Hence, computing Rp reduces to the solution of 
the simple equation: 

J2 



2E + 2 



Ap + 0.5M 



2Bp + — = 0. 



(C.l) 



R ' • i?2 

Needless to say, the computation of the apo-centre radius 
is very similar. 

Appendix D: Tests for spurious relaxation 

To measure the amplitude of the spurious relaxation due 
to the fact that we represent a smooth potential $ by a 
set of super-stars, i.e. of spherical shells with zero thick- 
ness, we carried out a few simulations in which relaxation 
was switched off. In such cases, the algorithm reduces to 
moving the super-stars on their orbits again and again. 
Ideally, the structure of the cluster should not display any 
evolution but statistical fluctuations. Figure. D.l shows 
the results of such tests realized with 4000 and 64000 
particles. The computations where stopped when the to- 
tal number of steps divided by the number of super-stars 
reached 4000. For comparison, our core-collapse si mula- 
tion for a Plummer model with 2 x 10^ particles (Sect. 7.2) 
required an average of 2300 steps per super-star. From 
these relaxation-free tests, we conclude that the amount 
of spurious relaxation is negligible if the number of super- 
stars is of order 16 000 or larger. 

However, given that most CPU time is spent in tree 
traversals and that node access probabilities (proportional to 
St^^) are highly larger in the centre than in the outer parts 
[Stmax / Stmin S> 100), it IS quite unfortunate that our "lazy" 
updating method keeps lengthening the search path to those 
most active central particles. For the same reason . build i ng the 



Wirth 



tree as an "optimal' binary search tree ( Knuth 197S ; ^ ^ 

1976|) instead of a balanced one would certainly turn out to be 
a valuable improvement. 
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Fig. C.l. Binary tree strueture con- 
taining 32 Super-stars (squares). 12 
Super-stars have been evolved, leav- 
ing empty nodes (hexagons) . In this 
diagram, the horizontal position of a 
node reflects the radius of the asso- 
ciated super-star. The logical struc- 
ture of this tree is stored in the ar- 
rays whose content is displayed on 
the right. 
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a) Search for a super-star 
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b) Insertion of a super-star 
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c) Removal of a super-star 
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d) Peri-centre determination 

Begin 




/jp (not found) 

*:visit <- *:Adam (root node) 

j ^ 1, ^ ^ 0.5M, B ^ 



K-isJt ^ 


8 + 


"''Visit 




^visit ^ 


-Aa 


H (SA-isit 




-f^visit ^ 


- B - 


H ^-Byisit — 





^Yes— (QCRv|t)>0 
I No 



- Yes — ^3§(-^™it) < )_ No - 



Candidate memorisation 

kp i ^visit: 'p ^ 'visit 

^r- Aisit, -Bp <— Svisit 



Is there a left son 



j Yes 


'i ' + 


(^^vislt 


A^^A 




^'visit 


^'left son 



No- 



[s there a right sor 




, Yes 


B 

^visit ^ 


1- *-Bvisit 
'i' right son 



kp > 



Yes 



,No 



Pericenter < Ri \/ node ? 
jp ■(- 0, ^p <- 0.5M 
Bp ^ i3 + iSrisit 



(^iiid) 

Fig. C.2. Flow charts for elementary binary tree routines, a) Search for super-star with rank i; on exit A and B are 
its potential coefficients (see Sect. 5.1). b) Insertion of a super-star of mass Mss at radius Rss- c) Extraction of the 
super-star attached to node k. d) Peri-centre determination for a super-star of mass AI. See text for the definition of 
Q{R). On exit, fcp is the number of the last node with radius smaller than the peri-centre radius (in order of increasing 
radius), ip the corresponding rank and ^p. Bp, the potential coefficients at peri-centre. 



26 



Nss = 4000 
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Nss = 64 000 
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Fig. D.l. Tests for spurious re- 
laxation in Plummer models with 
4000 (left) and 64000 (right) super- 
stars. Relaxation was switched off. 
The top panels show the evolution 
(or absence thereof) of Lagrangian 
radii. The bottom panels present the 
evolution of central properties; den- 
sity (p), potential (<&) and veloc- 
ity dispersion (ct). These quantities 
have been normalized by the aver- 
age of the 20 first values of each 
sequence and the curves have been 
smoothed for the sake of clarity. 
Note that for A^ss = 4000, some spu- 
rious evolution happens while it is 
nearly suppressed for Ngs = 64000. 
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